Forêts aléatoires : Interprétabilité avec
les packages randomForestExplainer et pdp
Ruth Jeyaranjan
2024-02-21
Introduction
Dans le monde de la data science et de l’IA, on parle souvent d’apprentissage automatique, de modèles de boîtes noires ou encore d’interprétabilité. A partir du 20e siècle, les chercheurs ont commencé à élaborer des méthodes permettant de comprendre le comportement et l’interprétation de ces algorithmes fonctionnant de manière opaque. Mais si ça fonctionne, pourquoi l’interprétabilité est-elle si importante ? Nous allons détailler cette notion à travers un exemple.
Supposons qu’une banque utilise un modèle d’apprentissage automatique pour décider d’accorder ou non un prêt à un individu. Ce modèle prend en compte des caractéristiques telles que l’âge, le nombre d’enfants, le statut matrimonial et d’autres variables pour estimer la probabilité qu’un individu rembourse correctement le prêt. Si cette probabilité est de 70 %, comment savoir si cette décision est juste ? Quels facteurs influencent ce pourcentage ? Le statut matrimonial ? Le nombre d’enfants ? Est-ce que certaines caractéristiques sont surévaluées en raison des tendances dans les données d’apprentissage ? Et qu’en est-il des individus dont le cas est atypique ? Peuvent-ils être injustement rejetés en raison de leur singularité ? L’aspect boite noire de l’algorithme soulève une multitude questions.
A travers ce document, nous allons aborder les techniques
d’interprétabilité des forêts aléatoires avec les packages
randomForestExplainer et pdp.
Les données
Les particules suspendus dans l’air proviennent de diverses origines, aussi bien naturel qu’anthropique. Ces particules ont des compositions chimiques très variable. L’observatoire de la qualité de l’air en Haute-Normandie, Air Normand, mesure la qualité de l’air tout les quart d’heure à travers son réseau d’une douzaine de stations. Ces données sont accesibles via le package VSURF.
Nous allons utiliser les données enregistrés par la station
gui.
Description du package VSURF : Le package VSURF est
basée sur les forêts aléatoires. + description
suppressMessages(library(VSURF))# chargement du jeu de données gui
data("gui", package = "VSURF")Le jeu de données est composé de 18 variables qui sont les suivantes :
PM10 : Concentration moyenne journlière de PM10, in μ g/m^3
NO, NO2, SO2 : Concentration moyenne journlière de NO, NO2 , SO2, in μ g/m^3
T.min, T.max, T.moy : Température minimale, maximale et moyenne journaliers, en degrés Celsius
DV.maxvv, DV.dom : Vitesse maximale journaliers et direction du vent dominante, en degrés (pour la direction du vent, 0 degré correspond au Nord)
VV.max, VV.moy : Maximum et moyenne journaliers de la vitesse du vent, en m/s
PL.som : Pluies journalières en mm
HR.min, HR.max, HR.moy : Minimum, maximum et moyenne journaliers de l’humidité relative, en %
PA.moy : Pression atmosphérique journalière, en hPa
GTrouen, GTlehavre : Gradient de température journalier, en Celsius
Ces données seront utilisés dans le cadre de la prévision de la
pollution de l’air par les particules en suspension PM10 en fonction des
autres variables météorologiques. La variable cible est
PM10.
| PM10 | NO | NO2 | SO2 | T.min | T.max | T.moy | DV.maxvv | DV.dom | VV.max |
|---|---|---|---|---|---|---|---|---|---|
| 13 | 42 | 48 | 8 | -0.6 | 4.4 | 1.173913 | 160 | 180.0 | 8 |
| 13 | 22 | 34 | 7 | -1.3 | 5.6 | 2.012500 | 20 | 22.5 | 7 |
| 21 | 44 | 55 | 13 | -4.2 | -0.7 | -2.933333 | 40 | 180.0 | 3 |
| DV.dom | VV.max | VV.moy | PL.som | HR.min | HR.max | HR.moy | PA.moy | GTrouen | GTlehavre |
|---|---|---|---|---|---|---|---|---|---|
| 180.0 | 8 | 5.347826 | 11 | 93 | 97 | 96.21739 | 1010.930 | 0.7 | -0.3 |
| 22.5 | 7 | 4.958333 | 3 | 78 | 96 | 87.87500 | 1016.975 | -0.5 | 0.0 |
| 180.0 | 3 | 2.250000 | 0 | 63 | 91 | 82.62500 | 1024.088 | 2.4 | 0.1 |
Traitement des données
summary(gui)gui <- subset(gui,!is.na(PM10))Toutes les observations où la variable cible est nulle ont étés supprimés car ces observations n’apporterons aucune informations. Il reste encore 49 observations avec des données manquantes, ce qui représente 4.46% des données d’origine.
Nous allons créer une sous version de gui sans données manquantes que l’on appelle gui_without_na.
gui_without_na <- na.omit(gui)Analyse rapide des données
Regardons la distribution de la variable cible :
Fig. 1 Distribution de la variable cible
Distribution des co-variables en fonction de PM10 :
par(mfrow = c(3, 3), mar = c(2, 2, 2, 0) + .1)
for (i in 2:18) {
plot(gui[,i] ~ gui$PM10, main = names(gui)[i])
}
Division des données en données d’apprentissage et de test
Description du package Caret: Le package caret
(abréviation de Classification And REgression Training) est un ensemble
de fonctions qui tentent de rationaliser le processus de création de
modèles prédictifs. Cette librairie sera utilisé dans le cadre de
l’échantillonage des données en échantillons train et test avec la
fonction createDataPartition().
suppressMessages(library(caret))Les données d’entraînement permettent d’entraîner le modèle
d’apprentissage, cela signifie que que le modèle va apprendre sur ces
données. Les données test permettent de tester la performance du modèle
avant de le déployer. On utilise souvent 70 à 80% pour entrainer le
modèle et 20 à 30% pour le tester. Si on utilise toutes les données
disponibles pour l’apprentissage du modèle, on ne pourra pas tester sa
performance. Ici, le ratio sera de 30/70.
set.seed(123)
partition <- createDataPartition(y=gui_without_na$PM10, p=0.7, list=FALSE)
training <-gui_without_na[partition,]
testing <- gui_without_na[-partition,]A. Construction d’une forêt aléatoire
RandomForest est un algorithme de classification composé de nombreux arbres de décision. Elle utilise l’échantillonnage et le caractère aléatoire des caractéristiques lors de la construction de chaque arbre individuel pour créer une forêt d’arbres non corrélés dont la prédiction par comité est plus précise que celle d’un arbre individuel.
Description du package randomForest : hehe
suppressMessages(library(randomForest))Forêt aléatoire simple
Nous allons commencer par construire une forêt aléatoire simple, sans toucher aux hyperparamètres et établir le RMSE (Root Mean Squared Error) de référence.
La formule pour calculer le RMSE est :
\[ \text{RMSE} = \sqrt{\frac{\sum_{i=1}^{n} (y_i - \hat{y}_i)^2}{n}} \]
set.seed(123)
default_RF <- randomForest(PM10 ~ NO+NO2+SO2+T.min+T.max+T.moy+DV.maxvv+DV.dom+VV.max+VV.moy+PL.som+HR.min+HR.max+HR.moy+PA.moy+GTrouen+GTlehavre, data = training)
print(paste("RMSE =",round(sqrt(mean((testing$PM10-predict(default_RF, testing))^2)),2)))## [1] "RMSE = 5.95"Nous obtenons un RMSE de référence de 5,95. Cela signifie qu’en moyenne la valeur prédite varie de plus ou moins 5,96 par rapport à la valeur réelle.
Nous allons aussi calculer la valeur du R² de référence.
La formule pour calculer R² est :
\[ R^2 = \left(\frac{\sum_{i=1}^{n} (y_i - \bar{y})^2 - \sum_{i=1}^{n} (y_i - \hat{y}_i)^2}{\sum_{i=1}^{n} (y_i - \bar{y})^2}\right) \]
r_squared=cor(testing$PM10, predict(default_RF, testing))^2
print(paste("R^2 =", round(r_squared, 3)))## [1] "R^2 = 0.644"Nous obtenons un R² de 0.64. Ce qui signifie que le taux de précision du modèle est de 64%.
Forêt aléatoire avec les paramètres ajustés
Les paramètres classiques du modèle sont appris au cours de l’aprentissage du modèle, mais les hyperparamètre doivent être définis auparavant. Dans le cas de RandomForest, les hyperparamètres pouvant être définis sont le nombre d’arbres de décision dans la forêt (ntree) et le nombre de caractéristiques prises en compte par chaque arbre lors de la division d’un nœud (mtry).
Le réglage de ces paramètres est plutôt expérimental que théorique, nous allons essayer de nombreuses combinaisons différentes et évaluer, comparer les résultats pour trouver les réglages optimaux. Dans cette partie, nous allons chercher la meilleure valeur de ntree avec do.trace et mtry en entraînant plusieurs modèles.
Choix du nombre d’arbres
Nous allons entrainer un modèle de forêt aléatoire pour prédire PM10
en fonction des autres variables avec l’argument do.trace
de randomforest. Le modèle construit 500 arbres de décision et affiche
des informations sur l’avancement de la construction de l’arbre toutes
les 25 itérations. Ensuite, le nombre d’arbres optimal sera choisi en
trouvant l’erreur OOB minimale.
set.seed(123)
RF_DoTrace <- randomForest(PM10~., data=training, ntree = 500, do.trace=25)## | Out-of-bag |
## Tree | MSE %Var(y) |
## 25 | 47.46 51.28 |
## 50 | 44.35 47.92 |
## 75 | 42.22 45.62 |
## 100 | 41.68 45.03 |
## 125 | 41.4 44.73 |
## 150 | 41.21 44.53 |
## 175 | 41.16 44.47 |
## 200 | 41.24 44.56 |
## 225 | 41.14 44.45 |
## 250 | 41.03 44.33 |
## 275 | 40.87 44.16 |
## 300 | 40.78 44.06 |
## 325 | 40.72 44.00 |
## 350 | 40.81 44.10 |
## 375 | 41.05 44.36 |
## 400 | 41.1 44.41 |
## 425 | 41.09 44.40 |
## 450 | 41.02 44.33 |
## 475 | 41.04 44.34 |
## 500 | 41.01 44.31 |On remarque que l’erreur diminue de 25 à 325 arbres puis augmente de 350 à 500 arbres. D’après cette sortie, le meilleur nombre d’arbre est de 325.
Choix du nombre de noeuds
Nous allons entrainter un modèle de forêt aléatoire (ntree=325) pour chaque valeur de mtry entre 1 et 17 et calculer l’erreur OOB. Puis, nous allons représenter les erreurs sous forme de nuage de points.
nbvars <- 1:(ncol(gui) - 1)
oobsMtry = c()
for (nbv in nbvars){
set.seed(123)
RF <- randomForest(PM10~., training, ntree = 325, mtry = nbv)
oobsMtry = c(oobsMtry,sqrt(mean((testing$PM10-predict(RF, testing))^2)))
}plot(nbvars, oobsMtry)
La valeur optimale de mtry est 12 car elle conduit à l’erreur OOB minimale.
On refait un modèle avec les bons hyperparamètres et on estime ses performances.
set.seed(123)
tuned_RF <- randomForest(PM10 ~ ., ntree = 325, mtry = 12, data = training, importance = TRUE)
print(paste("RMSE =",round(sqrt(mean((testing$PM10-predict(tuned_RF, testing))^2)),2)))## [1] "RMSE = 5.92"r_squared=cor(testing$PM10, predict(RF_DoTrace, testing))^2
print(paste("R^2 =", round(r_squared, 3)))## [1] "R^2 = 0.644"Conclusion :
En ajustant les hyperparamètres, nous obtenons un RMSE légèrement inférieur (5,92 < 5.96) mais un R^2 aussi légèrement inférieur (0,642 < 0.644). Les deux modèles montrent des performances très similaires, autant garder le modèle le plus simple.
Le modèle final qui sera utilisé dans la suite du document est le modèle sans fine tuning.
Importance des variables avec randomForest
varImpPlot(default_RF, sort =T, main = 'RF - Variable importance ')
Les variables les plus importantes sont NO2, NO, GTrouen, VV.moy, PA.moy, GTlehavre et SO2.
C. Interprétabilité avec le package randomForestExplainer
Description du package randomForestExplainer :
fonctions qui calculent un ensemble de mesures d’importance variable, y compris celles déjà utilisées et quelques nouvelles telles que le nombre de fois qu’un prédicteur a été utilisé pour diviser le nœud racine d’un arbre dans la forêt.
grande variété de possibilités de visualisation de la forêt qui aident à comprendre sa structure et à évaluer l’importance des variables.
fonction wrapper qui utilise toutes les capacités du package pour une forêt aléatoire donnée et produit un rapport HTML résumant les résultats.
une fonction permettant de sortir un fichier html avec tout les graphiques primaires du package (explain_forest()).
suppressMessages(library(randomForestExplainer))Minimal depth distribution
La fonction min_depth_distribution permet de calculer la
profondeur des prédicteurs pour chaque arbre et rassembler les résultats
pour n arbres dans un seul dataframe et calcule le minimum de la
profonduer de chaque variable dans chaque arbre.
#min_depth_frame = min_depth_distribution(default_RF)
#save(min_depth_frame, file = "min_depth_frame.rda")
load("min_depth_frame.rda")
knitr::kable(head(min_depth_frame, n=10), caption = "Distribution de minimal_depth")| tree | variable | minimal_depth |
|---|---|---|
| 1 | DV.dom | 2 |
| 1 | DV.maxvv | 6 |
| 1 | GTlehavre | 5 |
| 1 | GTrouen | 0 |
| 1 | HR.max | 5 |
| 1 | HR.min | 7 |
| 1 | HR.moy | 5 |
| 1 | NO | 1 |
| 1 | NO2 | 3 |
| 1 | PA.moy | 4 |
La fonction plot_min_depth_distribution utilise les
informations sur la profondeur minimale de chaque variable dans chaque
arbre de la forêt pour représenter graphiquement la distribution
discrète de la profondeur minimale.
plot_min_depth_distribution(min_depth_frame)
Les varibles incluses dans le graphiqud sont NO, NO2, GTrouen, PA.moy,VV.moy, GTlehavre, SO2, T.moy, T.max et T.min. Les variables NO et NO2 (1.95 et 1.97) ont des profondeurs mimales moyennes plus basses, ce qui signifie qu’elles sont plus influentes dans la prédiction du modèle. En revanche, es variables T.min, T.max et T.moy (3.67, 3.79 et 3.79) ont des profondeurs minimales moyennes plus élevées, ce qui signifie qu’elles sont moins influentes que les autres, dans ces 10 variables.
En somme, ce graphique permet de faire de l’analyse de données et de sélectionner des variables dans le cadre de la modélisation des données.
Variable importance
La fonction plot_multi_way_importance utilise les
données générés (datafrme contenate plusieurs mesures d’importance de
variables) par la fonction measure_importance et trace deux
ou trois mesures d’importances de variables les unes contre les
autres.
#importance_frame = measure_importance(default_RF)
#save(importance_frame, file = "importance_frame")
load("importance_frame")plot_multi_way_importance(importance_frame, size_measure = "no_of_nodes")
Les points sont allignés le long d’une diagonale, ce qui suggère une corrélation positive entre le nombre de noeuds et la profondeur minimale moyenne. On identifie les variables les plus importantes dans la zone supérieur droite du graphique : NO et NO2. Ce qui correspond bien avec graphique “minima depth distribution”.
plot_multi_way_importance(importance_frame, size_measure = "p_value")## Warning: Using alpha for a discrete variable is not advised.
La fonction plot_importance_ggpairs utilise les données
générés (dataframe contenant plusieurs mesures d’importance de
variables) par la fonction measure_importance et représente
graphiquement les corrélations entre les différentes mesures
d’importance des variables.
plot_importance_ggpairs(importance_frame)
Les mesures de node_putity_increase et mean_min_depth sont fortement corrélés (-0.933). La forte corrélation négative entre ces deux variables suggère que les nœuds améliorant la pureté sont situés à des profondeurs plus élevées dans les arbres de décision, alors que les variables moins profondes ont une influence moindre sur la pureté des nœuds.
Variable interactions
La fonction plot_min_depth_interactions prend en entrée
le résultat de la fonction min_depth_interactions qui
contient des informations sur la profondeur minimale moyenne
conditionnelle pour chaque interaction entre variables, calculé de
différentes manières. Elle trace ensuite la profondeur minimale moyenne
conditionnelle pour un nombre donné des interactions les plus
fréquentes, ainsi que la profondeur minimale moyenne de la deuxième
variable dans chaque interaction afin de comparer. Cette représentation
graphique permet donc d’identifier les relations les plus importantes et
de comprendre comment les différentes variables interagissent pour
influencer les prédictions du modèle.
#(vars <- important_variables(importance_frame, k = 5, measures = c("mean_min_depth", "no_of_trees")))
#interactions_frame <- min_depth_interactions(default_RF, vars)
#save(interactions_frame, file = "interactions_frame.rda")
load("interactions_frame.rda")
knitr::kable(head(interactions_frame[order(interactions_frame$occurrences, decreasing = TRUE), ]), caption = "Visualisation de intercations_frame généré par la fonction min_depth_interactions")| variable | root_variable | mean_min_depth | occurrences | interaction | uncond_mean_min_depth | |
|---|---|---|---|---|---|---|
| 48 | PA.moy | NO2 | 1.863839 | 448 | NO2:PA.moy | 3.118 |
| 13 | GTlehavre | NO2 | 2.407741 | 446 | NO2:GTlehavre | 3.278 |
| 38 | NO | NO2 | 2.106402 | 446 | NO2:NO | 1.954 |
| 73 | T.moy | NO2 | 2.068085 | 445 | NO2:T.moy | 3.672 |
| 42 | NO2 | NO | 1.951982 | 444 | NO:NO2 | 1.974 |
| 43 | NO2 | NO2 | 2.105290 | 443 | NO2:NO2 | 1.974 |
plot_min_depth_interactions(interactions_frame)
Les cinq interactions les plus fréquentes sont NO2:PA.moy, NO2:GTlehavre, NO2:NO, NO2:T.moy et NO:NO2.
Les quatres interactions les plus importantes sont NO2:HR.max, NO:HR.max, NO:DV.dom et NO2:DV.dom.
Seule l’interaction NO2:PA.moy se situe au même niveau que la ligne rouge, cette interaction a une profondeur minimale moyenne conditionnelle similaire à la moyenne des interactions. Toutes les autres interactions ont une profondeur minimale moyenne conditionnelle plus élevée.
Prédiction de la forêt sur une grille
Nous avons calculer les interactions, maintenant nous pouvons nous intéresser à la manière dont la prédiction de la forêt dépend des valeurs des deux variables qui constituent l’interaction.
La fonction plot_predict_interaction permet de tracer la
prédiction de notre forêt sur une grille de valeurs pour les composants
de chaque interaction. Nous allons appliquer cette fonction aux quatres
interactions les plus importantes.
NO2:HR.max
plot_predict_interaction(default_RF, gui_without_na, "NO2", "HR.max")
Le graphique montre une teinte bleue prononcée à gauche, et des nuances plus claires et du rouge à droite. Des valeurs plus élevés de NO2, conduisent à des prédictions plus élevés de PM10.
NO2:HR.max
plot_predict_interaction(default_RF, gui_without_na, "NO", "HR.max")
Le graphique montre une teinte bleue prononcée à gauche, et des nuances rouges à droite. Des valeurs plus élevés de NO, conduisent à des prédictions plus élevés de PM10.
D. Interprétabilité avec le package pdp
Description du package pdp : Ce package a pour but de
construire des courbes de dépendance partielle (c’est-à-dire des effets
marginaux) à partir de différents types de modèles d’apprentissage
automatique en utilisant R.
library(pdp)Les courbes de dépendance partielle (PDP) ont en effet été introduits par Friedman en 2001, principalement dans le but d’interpréter les modèles d’apprentissage automatique complexes, en particulier ceux qui dépassent la simple régression linéaire. Le graphique analyse comment la variation du prédicteur i influence les prédictions du modèle en moyennant les autres prédicteurs. En examinant comment les prédictions du modèle changent à mesure que les valeurs de i varient, on peut comprendre l’impact de cette variable sur les résultats du modèle, facilitant ainsi l’interprétation de son importance dans la prédiction.
Nous allons construire des graphiques partiels pour les variables explicatives qui sont les plus importantes dans notre modèle. Chaque graphique illustre la relation entre une variable explicative et la prédiction du modèle pour PM10, en tenant compte des autres variables dans le modèle.
partial(default_RF, pred.var = "NO", plot = TRUE)
partial(default_RF, pred.var = "NO2", plot = TRUE)
partial(default_RF, pred.var = "SO2", plot = TRUE)
partial(default_RF, pred.var = "T.min", plot = TRUE)
partial(default_RF, pred.var = "T.max", plot = TRUE)
partial(default_RF, pred.var = "T.moy", plot = TRUE)
partial(default_RF, pred.var = "VV.moy", plot = TRUE)
partial(default_RF, pred.var = "PA.moy", plot = TRUE)
partial(default_RF, pred.var = "GTrouen", plot = TRUE)
partial(default_RF, pred.var = "GTlehavre", plot = TRUE)
La valeur prédite de PM10 augmente quand les NO et NO2 augmentent. Elle augmente aussi quand SO2 augmente. La valeur prédite augmente quand T.min diminue en dessous de -5 ou augmente au dessus de 10. Elle augmente quand T.max va en dessous de 2 et augmente entre 13 à 30. Elle augmente quand T.moy est en dessous de 0 et augmente entre 10 et 25. Elle décroit strictement qd vv.moy est entre 1 et 2, puis stagne entre 2 et 10 à une valeur basse.
Les variables NO et NO2 ont un effet positif sur la valeur prédite de PM10, ce qui signifie que des niveaux plus élevés de ces polluants sont associés à des valeurs prédites plus élevées de PM10. De même, l’augmentation de SO2 est également associée à une augmentation de la valeur prédite de PM10.
Pour les variables liées à la température (T.min, T.max, T.moy), des tendances complexes sont observées : la valeur prédite de PM10 augmente lorsque T.min est inférieur à -5 ou supérieur à 10, diminue lorsque T.max est inférieur à 2, puis augmente entre 13 et 30. De plus, la valeur prédite augmente lorsque T.moy est inférieur à 0 et augmente entre 10 et 25. Quand il fait très froid ou très chaud, le taux de PM10 prédit augmente.
En ce qui concerne la variable vv.moy, la valeur prédite de PM10 décroît strictement lorsque vv.moy est entre 1 et 2, puis stagne entre 2 et 10 à une valeur basse.
Il est aussi possible de faire une représentation tridimensionnelle
qui montre la relation conjointe de deux variables explicatives avec la
prédiction du modèle, grâce à plotPartial().
pd <- partial(default_RF, pred.var = c("NO", "NO2"), chull = TRUE)
plotPartial(pd, main = "représentation tridimensionnelle entre NO et NO2 avec la prédiction")
Conclusion
En utilisant le package randomForestExplainer, nous
avons exploré la distribution des profondeurs minimales des variables et
évalué leur importance à l’aide de divers graphiques. Les variables les
plus influentes identifiées incluent la Concentration moyenne
journalière de PM10, NO, NO2, SO2, la Température minimale, maximale et
moyenne journalières, la Vitesse maximale et moyenne du vent, les Pluies
journalières, l’Humidité relative minimale, maximale et moyenne, la
Pression atmosphérique journalière, et le Gradient de température
journalier pour les stations GTrouen et GTlehavre. En examinant les
interactions entre les variables avec plot_min_depth_interactions, nous
avons visualisé comment elles affectent les prédictions du modèle.
Enfin, en utilisant plot_predict_interaction, nous avons exploré
l’impact des niveaux de polluants sur les prédictions de qualité de
l’air.
En utilisant le package pdp, nous avons utilisé des
courbes de dépendance partielle (PDP) pour interpréter un modèle de
forêt aléatoire appliqué à des données environnementales. Les PDP
permettent de visualiser comment la variation d’une variable explicative
influence les prédictions du modèle tout en tenant compte des autres
variables. Nous avons construit des graphiques partiels pour les
variables les plus importantes telles que la Concentration moyenne
journalière de NO, NO2, SO2, les Températures minimale, maximale et
moyenne journalières, la Vitesse maximale et moyenne du vent, la
Pression atmosphérique journalière, et le Gradient de température
journalier pour les stations GTrouen et GTlehavre.